Modelling a Dune Field 



A. R. Lima 1 , G. Sauermann 1 ' 2 , H. J. Herrmann 1 ' 2 and K. Kroy 1 

1) PMMH, Ecole Superieure de Physique et Chimie Industrielles (ESPCI), 10, rue 
{NJ • Vauquelin, 75231 Paris, Cedex 05, France 

2) ICA-1, University of Stuttgart, Pfaffenwaldring 27, 70569 Stuttgart, Germany 

(N ' 

Abstract 

00 

We present a model to describe the collective motion of barchan dunes in a field. 
Our model is able to reproduce the observation that a typical dune stays confined 
within a stripe. We also obtain some of the pattern structures which ressemble those 
observed from aerial photos which we do analyse and compare with the specific field 
■ of Laayounne. 
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1 Introduction 

a 

A wide variety of dune shapes can be found in deserts, in coastal areas, on the 
sea-bottom [1], and even on Mars [2] and over one hundred different types of 
dunes have been classified by geomorphologists and their ancient predecessors 
so that some of them bear old arabic names. The shapes depend mainly on 
the amount of available sand and on the change in the direction of the wind 
over the year. If the wind blows steadily from the same direction throughout 
the year and there is not enough sand to cover the entire area, dunes called 
barchans develop having the shape of a crescent. They typically appear in 
huge fields of several hundred kilometers length stretching along coastlines 
and driven by trade-winds. A picture from a section of such a field is shown 
in Fig. 1. 

Measurements of the relationship between height, width and length have been 
performed for instance in Peru or Morocco [3-8] and led to a good under- 
standing of the barchan shape. In particular we know that the barchans move 
without changing their shape with a velocity which is inversely proportional to 
their height. However, the dynamics of dunes in a dune field has not yet been 
studied very much. In particular, there seems to be a contradiction between 
the fact that the barchans are not all of equal height and the stability of the 
dune field since when the dunes have different velocities they will eventually 
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Fig. 1. Aerial photograph of a barchan dune field close to Laayounne in Morocco. 
The photo clearly shows a correlation between the dunes in the lateral direction. 

bump into one another. We do not yet understand the underlying stabilizing 
mechanism. Looking at the images in more detail also other questions arise. 
Are the relative positions of the dunes correlated? Do they form patterns? 
When dunes are getting closer they can interact through the non-local wind 
field or the sand flux between the dunes. In this work, we propose a simple 
model that defines the interaction between the dunes caused by the inter-dune 
sand flux. We neglect further effects of the wind field. The model predicts the 
time evolution of the entire field. It can reproduce the spatial correlations be- 
tween the dunes that have been observed in the field, cf. Fig. 1. Furthermore, 
the model shows that the entire dune field does not broaden but stays focused. 
However, the stability of the dune field depends on the sand flux balance of 
a single dune, which is not well known so that we will have to make some 
simplifying assumptions. 



2 Wind and sand flux 

Sand is transported if the wind speed exceeds a certain threshold velocity. 
In equilibrium the amount of transported sand or the flux of sand q(u*) can 
be expressed as a function of the shear velocity u* [9]. The shear velocity 
u* characterizes the turbulent boundary layer of the atmosphere and thus its 
well known logarithmic velocity profile v(z) = u*/ /tin z/z , where v(z) is the 
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wind speed at height z over the ground, n ~ 0.4 the universal von Karman 
constant for turbulent flow, and z the roughness length of the surface. For 
shear velocities well above the threshold, e.g. = 0.5 ms 4 , the simplest sand 
flux relation, proposed by Bagnold [10], should give a reasonable prediction 
for the sand flux q, 

q = C B ^ul (1) 
9 

where Cb ~ 2 is a phenomenological constant [10-12], p a = 1.225 kg m -3 the 
density of air, and g = 9.81ms -1 the acceleration due to gravity. Using the 
bulk density p s = 1650 kg m -3 of a typical sand we can rewrite Eq. (1) into a 
relation for the volume flux $, 

$ = (2) 

Ps 

For the discussion of the dune field dynamics we chose the shear velocity 
over the flat ground m*o = 0.5 ms -1 and obtain a volume sand flux $ °f 
0.031 m 2 s -1 . Furthermore, we assume that the sand velocity v s is proportional 
to the shear velocity u* and choose a typical value of v s = 5 ms -1 . 



3 Barchan dunes 



A dune is characterized by a single scalar quantity, e.g. its width w, and its 
position (x, y) in the dune field. Further properties such as the height h, 

h = aw (3) 

or volume V, 

V = bw 3 (4) 

can be calculated using simple phenomenological formulas. Here, the param- 
eters a = 0.09 and b = 0.05 have been obtained from field measurements 
[8]. 

The migration velocity of a dune depends on its height h and the sand 
flux $ b over the brink of the dune [10]. The latter can be calculated through 
Eq. (1) if the shear velocity at the brink is known. However, the dune itself 
disturbs the wind field such that it becomes stronger above the dune than 
over a flat ground. The ratio between the shear velocity over a flat ground u*o 
and the shear velocity at the brink of the dune is called the "speed-up" n. 
This speed-up depends on the shape of the dune and mainly on the average 
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Fig. 2. Sketch of a barchan dune. 



windward slope. A typical value for the speed-up over a sand dune is 1.4 
[13,14]. Finally, we can write for the dune velocity Vd, 

$6 <5>on 

v d = -r = • 5 

h aw 

We want to note that Eqs. (3) and (4) implicitly assume scale invariant dune 
shapes, which is an approximation that is only valid for large dunes, much 
larger than the minimal dune size [15,16]. 

The growth and shrinkage of a dune depends on the balance of in- and out- 
going sand fluxes. The sand flux that arrives from a upwind location, for 
instance a beach or other dunes, feeds the dune with new sand, whereas at 
the tip of the horns sand leaves the dune. The change of volume Vi due to the 
in-flux $(y) can be calculated by integrating over the dunes width, 

Vi= / $(y)dy. (6) 
Jo 

The sand that leaves the dune at the horns gives rise to a loss of volume V D 
that is modeled by the following expression, 

V = w %(a + axw) , (7) 

where ao and a\ are parameters that depend on the complex processes oc- 
curring on the dune's surface and are not known a priori. Recent numerical 
calculations of the sand outflux from barchan dunes suggest that the outflux 
is nearly proportional to the dunes width and motivated the expression above 
[16]. 

For a stable dune, the gain and loss of volume must exactly compensate Vi — 
V Q = 0. The average inter-dune flux or influx ($) = / $ causes an average 
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Fig. 3. Sketch of a dune field. 



increase in volume, 

[V t ) = wf%, (8) 



where / defines the saturation of the inter-dune flux. Eqs. (7) and (8) define 
finally an average dune width (w), 

<«> = (9) 

Using the average dune width (w) and 7 = a x (w) //as new parameters, 
instead of ao and a±, we can rewrite Eq. (7), 

K = tt ./*„(i +7 ^M). (io) 

Finally, the change in volume V of a dune can be obtained from Eqs. (6) and 
(10), 

V = £^y)dy-wf%[l + 1 ^^^ (11) 

Typical values for / are of the order of 0.1, i.e. the inter-dune flux is far from 
being saturated. 



4 The barchan field 



A barchan field of length If and width Wf as depicted in Figure 3 is composed 
of many dunes that migrate in the x-direction, the direction of the wind. The 
center of mass (x,, yi) of a dune is used to localize the dune i in the field. 
The dynamics of the dune field is described by the temporal change of the 
variables Xi, and W{ which will be defined in the following. 



From Eqs. (11) and (4) we obtain the change of the width of a dune, which 
depends on the inter-dune flux, 



w 



ryi+Wi/2 ( w _ ( w \ ' 

/ *(y)dy- Wi f* + 

Jyi-Wi/2 \ {W) 



1/3 



(12) 



The velocity of the dune in the wind direction is given in Eq. (5) and depends 
only on the width of the dune, 



Xi = V d (Wi) = 



avjj 



(13) 



An inter-dune flux that is asymmetric in the lateral direction with respect to 
the center of mass yi of the dune gives rise to a change in lateral position y^ 



VWi/2 



where V is the volume of the dune i. 



(14) 



Dunes can interact either due to the inter-dune flux or due to coalescence. 
The first case is defined by the equations above and the inter-dune sand flux 
$(x, y) defined in the next section. For the latter case we define a characteristic 
distance of coalescence d co \- If two dunes become closer than this distance they 
are merged to a single dune. The mass is conserved by this process and the new 
position is calculated according to their centers of mass and relative volumes. 



4-1 The inter-dune sand flux in the field 



The flux of sand in the field, i.e. between the dunes, is a complex problem 
and causes an interaction between the dunes. Dunes can act here as sources 
and sinks for the sand at the windward foot and at the horns, respectively. 
Furthermore, the position and width of the dunes changes in time and thus 
their influence on the inter-dune sand flux. We propose the following model 
for the inter-dune flux: First, the sand is convected with the constant velocity 
v s in the wind direction. Secondly, a lateral movement is caused by diffusion 
and obeys a diffusion equation, 

M- D W = °- (15) 

where the time scale of the diffusion process is connected to the transport in 
wind direction by the sand velocity v s . Hence, we start solving the diffusion 
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equation at the influx-border x = of the field, where the sand flux is spatially 
constant. When integrating the solution forward in time by displacing its re- 
position we obtain the stationary solution of the two-dimensional inter-dune 
flux field 

^~ D w=°- (16) 



For each dune in the field, the flux $j is set to zero along the entire width of 
the dune (from y^ — Wi/2 to yi + u>j/2). Then, at the horns (i.w. the positions 
Hi — Wi/2 and yi + Wi/2) the outflux is imposed as two 5-functions. Finally, 
the new dune mass is calculated by the in- and out-flux balance, Eq. (12). 



4-2 Initial configuration and Boundary conditions 



The field is initialized with dunes at random positions until a certain density 
Pd is reached, which is a parameter of the model. The dune width is chosen 
between a minimum and a maximum width, w m i n < w < w max , having an 
average width (w) = (w min + w max )/2. From the initial density pd we can 
define the rate r with which new dunes enter the field at the in-flux boundary, 

r = Wp d v((w)), (17) 

where v((w)) is the dune velocity according to Eq. (5) and W is the width of 
the dune field, which is chosen to be only one third of the total simulation area, 
cf. Fig 3 where Wf = 3W. The width of the dunes that enter with the rate r 
at the influx boundary is chosen randomly as in the initial configuration. If a 
dune crosses the outflow boundary (at the right side of the field) it is removed 
from the simulation. 

The inter-dune flux of sand <3> is chosen to be constant between Wf/3 < y < 
2wf/3 and zero outside. The constant influx is chosen to be $ = /$ and thus 
corresponds to an inter-dune flux of a dune field with average dune width (w), 
cf. Eq. (8). In the direction perpendicular to the wind (^/-direction) periodic 
boundary conditions are used. 



5 Numerical Scheme 



The numerical scheme can be divided into several steps. First, the dune field 
is randomly initialized with dunes as discussed above. Then, an iterative cal- 
culation, forward in time, of the following steps is performed: 
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• New dunes are created with rate r at the left boundary, in the region Wf/3< 
y<2w f /3. 

• The new x-position of the dunes is calculated according to Eq. (13). 

• If two dunes are closer than the coalescence distance d co \ = 50 m they are 
merged. 

• The inter-dune sand flux is calculated from left to right on a discrete one di- 
mensional grid in ^/-direction according to Eq. (16). The iteration is started 
at x = with a constant value of $(y) = /$ f° r w f/3 < V < 2t«//3 and 
zero outside. When the calculation forward in the wind direction reaches a 
dune, the flux is set to zero along the dune width and two discrete 5-peaks 
are set at the horns of the dune as described above. Furthermore, the change 
in volume of the dune and its lateral movement according to Eqs. (12) and 
(14) are calculated. 

The time step between the iteration cycles above has been chosen to be AT = 
1/12 year. 



6 Results 

We are mainly interested in measuring how the dunes are distributed in the 
field, possible correlation between the positions of the dunes and understand- 
ing the confinement of the field on a stripe. We did simulations of our model 
for a field having a length of 50 km and a width of 12 km. The dunes enter the 
field only in the region between km 4 and km 8. The initial density of dunes 
in the field was chosen to be pf = 10 dunes/km 2 . We set the saturation of the 
flux to be / = 0.15 which corresponds to a flux of = 0.0046 m 2 /sec. 
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Fig. 4. Snapshot of the simulation of the field. 

In fig. 4 we see a typical simulation of the dune field. It stays confined inside 
a stripe. The dune density slowly decreases in direction of the wind. This 
can be seen quantitatively in Fig. 5 where the density averaged over the y 
direction is plotted over the x-axis (x and y directions are defined in fig. 4) 
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in a semilogarithmic plot. We recognize for large x a logarithmic decay of the 
form p oc In x. 

The situation seen in fig. 4 is a stationary state. As seen in fig. 6 the average 
size of the dunes fluctuates around a value that is constant in time and except 
for a short transient at early times also the number of dunes is essentially 
constant in time. 
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Fig. 6. Number and average size of dunes in the field 
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6. 1 The confinement of the field 



It is not clear up to now why real dune fields are typically confined on a stripe 
all the time. In our simple model, despite the fact that the diffusion of the 
sand tends to spread the sand to the borders, this behavior appears by itself. 
As seen in fig. 7, the sand flux decreases outside the central region of the field. 
Hence, dunes at the borders of that region receive less sand. The ones close to 
the center of the field receive more sand. This flux gradient acts like a "force" 
which pulls the dunes into the region where there is more flux, i.e. the central 
region where the other dunes are. 

A way to see this effect is the distribution of dunes and their mean size at the 
end of the field. This is shown in figure (8). In the inset we plotted the number 
of dunes that leave the field at the right boundary. Near to the center of the 
field the dunes have approximately the same size, however the ones near to 
the borders are smaller, since they receive less sand. Since the dunes in the 
borders are smaller, they are faster and this explains why we have more dunes 
leaving the field close to the borders than in the center. 
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Fig. 7. Averaged flux at the end of the field. 



6.2 Spatial correlations 

Looking at fig. 1 we can observe that the dunes tend to form lines, the horn 
of one dune pointing to the center of the back of another one. These patterns 
are apparently due to the coupling through the inter-dune sand flux and we 
will show how they are reproduced by our model. In order to quantitatively 
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Fig. 8. Average size and number of dunes (inset) which leave a field of 50 km. We 
can see that the field stays confined and that the dunes outside of the central regions 
tend to shrink. 



measure these structures, we will introduce two types of correlation functions 
which we will monitor as well in our model as also directly from a photograph 
of a real dune field in order to compare the patterns formed by our model to 
the observed ones. 

Let us define the longitudinal correlation as 



This function measures the correlation of dune positions in the wind direction. 

In fig. 9 we present the results for our model. To get better statistics we take 
the average over several thousand time-steps. We observe very pronounced 
peaks at regular distances of about 50 m showing a characteristic distance 
between the dunes in wind direction. 

It is also interesting to calculate the angular correlations in the positions of 
the dunes which can be quantified by the angular correlation function 
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Fig. 9. Correlations in the wind direction obtained from our model after averaging 
over time. 

where i and j are only chosen within the central part of the stripe (about 
one third of all dunes) in order to avoid boundary effects and for which the 
distance (xj — Xj) 2 + {Ui + Uj) 2 < (5km) 2 . In fig. 10 we see our results for C p (9) 
as obtained after averaging over all dunes in the central part of the field and 
several thousand time-steps. 
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Fig. 10. Polar plot of the angular correlation C p {9) between the dunes. 

We see that while most directions have no particular structure in the forward 
direction there are two favorite directions symmetrical around the exact wind 
direction (x-axis) of high correlation. In the wind direction itself the correlation 
is particularly low. The selected direction corresponds to the arrangements 
where the dune in front is just shifted by half a dune width, so that the horns 
directly point to the center of the dune in front. Such configurations can be 
recognized with the naked eye on the photo of fig. I. 
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7 Results from the experiments 



We analysed aerial photos of the dune field of Laayounne in Southern Morocco 
and obtained the position of 734 dunes as shown in fig. 11. In figs. 12 and 
13 we show our results for both correlation functions corresponding to the 
correlation functions obtained for the model in figs. 9 and 10. Comparing 
figs. 9 and 12 one sees that the measured longitudinal correlation does not 
show much periodicity or characteristic distance. We believe that this lack of 
structure is rather due to insufficient statistics. After all the model gives the 
possibility to the average over time while the measurements are obtained from 
a single photo. 
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Fig. 11. Position of the dunes in the Barchan field of Laayounne obtained from an 
aerial photo. 
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Fig. 12. Correlations in the wind direction obtained from the positions of the dunes 
in fig. 11. 



Concerning the polar correlations, fig. 13 clearly shows the same two peaks 
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Fig. 13. Polar plot of the angular correlations C p (9) of the dune field of Laayounne. 

pointing in both directions with an angle of 9 from the wind direction. But 
while in the model 9 was about 8°, it is about 15° for the real data. 



8 Conclusions 

In this paper we introduced a simple model for a field of Barchan dunes. The 
dunes are characterized by their width and position in the field. The model 
includes a diffusion of the sand. Assuming a simple law for the inter-dune sand 
flux we could show that the simulated field is stably confined and presents the 
same kind of spatial correlations as real fields. 

An extension of this work could be to introduce variations in the wind direc- 
tion. This could help to explain more complex patterns observed in deserts. 
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